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ABSTRACT 

A nonlinear elastic force— displacement relationship is used to calculate the 
transient impact force and local deformation at the point of contact between impactor nnd 
target. The nonlinear analysis and transfer function capabilities of NASTRAN are used to 
define a finite element model that behaves globally linearly elastic, and locally nonlinear 
elastic to model the local contact behavior. 

Results are presented for two different structures: a uniform cylindrical rod 
impacted longitudinally; and an orthotropic plate impacted transversely. Calculated 
impact force and transient structural response of the targets are shown to compare well 
with results measured in experimental tests. 


INTRODUCTION 

Aerospace structures are subjected to impact loading from a variety of sources, 
including dropped tools, runway debris, and munitions. In some advanced materials, even 
low velocity impact can cause significant structural damage. Therefore, development of 
accurate means of calculating structural response due to impact loading can be of critical 
importance. In this paper, a computational technique is developed to predict the dynamic 
response of a structure to low velocity elastic impact. 

Structural damage due to impact invariably initiates in the immediate vicinity 
of the impact. Therefore, it is important that the local stress field in the region of contact 
be calculated accurately. Hertz [1] derived an elasticity— based force— displacement 
relationship that describes contact between two elastic bodies. The Hertzian contact law is 
given by: 


F = Kq" 


( 1 ) 


where 


F = elastic contact force 
K = contact stiffness 
n = exponent 
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and 


a — relative displacement (indentation) between impactor and target 
= Ui — ut (i = impactor, t = target) 

The exponent n was shown in reference [1] to have the value of 3 / 2 . In dynamic 
applications such as this; F, u, and a are all time— varying. 

During low velocity impact, where impact damage is confined to the area 
immediately around the point of contact, areas of the structure remote from the impact 
may still deform in a linear elastic manner. An efficient finite element model, therefore, 
would combine a linear elastic model of the global structure with a non— linearly elastic 
behavior at the point of contact. The nonlinear force— displacement relationship in 
equation (1) is incorporated into a linear elastic finite element model (MSC/NASTRAN 
transient solution 27, COSMIC/NASTRAN transient solution 9) by using a NASTRAN 
transfer function definition and nonlinear analysis capability. In the following section, the 
Hertz contact law is discussed in addition to a method of incorporating it into NASTRAN. 
Impact loading of two different structures is then analyzed. The first problem is a 
one-dimensional rod of uniform cross section impacted longitudinally. The second is an 
orthotropic plate under transverse impact. 


CONTACT LAW 

In reference [2] Hertz derived the force— displacement relationship for two 
spherical isotropic elastic bodies of radius rj, and r 2 in contact: 

F = K a 3 / 2 (2) 


where 


K 


4 

~ "T~ 


ri r 2 k t k 2 
r i + r 2 ki + k 2 


( 3 ) 


is the contact stiffness and 


k,. = i = 1,2 (4) 

1 — V • 
t 

where E^- and are the Young’s modulus and poisson’s ratio, respectively, and the 

subscripts 1 and 2 refer to each of the spheres. When a spherical impactor contacts a flat 
target, (3) simplifies to 


K 



ki k t 


( 5 ) 


where i and t represent the impactor and target respectively and the k t and ki are given by: 
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( 6 ) 



l — Vi 


In equation (2), a is the local indentation at the contact point, shown 
schematically in figure 1. We have: 

a = Ui-u t ( 8 ) 

where a is the relative local displacement between impactor and target at the point of 
contact. 

NASTRAN Implementation 

The non-linear local behavior was incorporated into the NASTRAN finite 
element model as follows: 

The impactor is modeled as a lumped mass just touching the target at t=0 and 
with an initial velocity towards the target. The difference between the displacement of this 
lumped mass and the displacement of the target is the indentation, or. The modeling of the 
contact between impactor and target is performed by utilizing the transfer function card, 
TF, and the nonlinear force card, NOLIN3. The TF card acts as a dynamic multipoint 
constraint, relating the displacement, velocity and acceleration of several independent 
degrees of freedom to a dependant degree of freedom. In the work discussed here, only 
displacement relationships were used. On the TF card coefficients of the following 
equation are specified [3]. 

n 

(B 0 + B lP + B 2 p2)u dep + £(A J 0 + A J lP + A J 2 p2)u J ind = 0 (9) 

j*i 

where 

B 0 , Bj, B 2 = the coefficients for the dependant degree of freedom 

A J 0 , A j, A 2 = the coefficients for the independent degrees of freedom 
u dep = ^he displacement of the dependant degree of freedom 

uj nd = the displacements of the independent degrees of freedom 
n = the number of independent degrees of freedom 
p = the differential operator ^ , and p 2 = — 

For this analysis, the equation would appear: 
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0 


( 10 ) 


(l.O)llextra point + l-0)Uimpactor + (l.O)utarget 

that is 

n = 2 

Bp B 2 , Aj, Aj = 0.0 (j = 1, n) 

B 0 = 1.0 

Aj = -10 

Aq = 1.0 

The resulting equation defines the indentation at every time step and assigns the value to 
an EPOINT. The EPOINT, or extra point, is used as a nonstructural variable that 
provides a place to store the value of the indentation. The EPOINT is provided as input 
to the NOLIN3 card. 

The NOLIN3 card is the means of applying the time— dependent nonlinear load 
based on the indentation. The NOLIN3 card has the form: 

s ( x ( 0 ) A » x (0 > 0 

P(t) = ■ (11) 

0 , x(t) < 0 

where 

P(t) = is the resulting nonlinear force 
S = is a scale factor 

x(t) = is the displacement or velocity of a degree of freedom 
A = is an amplification factor 

In modeling of the impact, we define x(t) to be the displacement of the EPOINT, S to be 
the Hertzian stiffness, and A to be 3/ 2 , a s given in equation (2). Recall that the 
displacement of the EPOINT is really the indentation as obtained from the TF card. The 
resulting function then has the form: 

K( Q (t)) 3 /’, o(t) > 0 

p (0 = ■ ( 12 ) 
0 , a(t) < 0 
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Note that when a is less than or equal to zero (ie. the target and the impactor are out of 
contact) then the force is also zero. Two NOLIN3 cards are used, one to apply the impact 
force to the target and the other to apply the same force to the impactor in the opposite 
direction of its initial velocity. This methodology allows the impactor to slow with 
increasing impact force and eventually to unload the target as the impactor begins to travel 
in the opposite direction, away from the target. 


RESULTS 


One Dimensional Rod 


The first problem analyzed is the longitudinal impact of a steel ball on a long 
aluminum rod of constant cross section. Geometry and material properties of the impactor 
and target are given in figure 1 . The problem was modeled using 144 1 — D rod elements 
with each grid point having a single longitudinal degree of freedom. Two more degrees of 
freedom were used to model the impact dynamics, resulting in a total of 147 degrees of 
freedom. A single lumped mass with an initial velocity was used to represent the impactor. 
The Hertzian force-displacement relationship in equation (1) was prescribed using the 
NASTRAN NOLIN3 card, as shown in the example input deck in the appendix. 

The calculated impact force history compares well with experimentally 
determined values [4], as shown in figure 2. The calculated strain response at the midpoint 
of the target bar is compared with measured values in figure 3. The sign reversal of the 
second pulse is caused by the reflected tensile stress wave generated by the incident 
compressive wave reaching the free end of the bar [5]. 

Some insight into the timing and the location of the impact— induced structural 
failure ran be gained by tracking the distribution of energy in the impactor and the target, 
as shown in figure 4. The energy balance can be expressed as: 

U t = KEj + SEi 4- KE t + SE t (13) 


where 


Ut = total energy in system 

2 

KEi = impactor kinetic energy = 72 mvj 

SEj = impactor strain energy = J F(a) da = 2 /s Ka 
KE t = target kinetic energy 

2 ( n _ num ber of elements) 

SEt = strain energy of target 


n - 1 


-X 


~Y~ m 


n - 1 




k i ( 


(n = number of elements) 


(14) 

(15) 

(16) 


(17) 
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The total energy in the system, Ut, is divided between the kinetic energy and 
the strain energy of the target and impactor in a time— varying manner. Because damping 
effects are not considered, the total system energy is constant and equal to the initial 
kinetic energy of the impactor. The strain energy of the impactor is non— zero only during 
the contact interval (0 < tC/L < 0.4, where t = time, L = the length of the bar, and C = 
the wave speed in the bar) and peaks when the contact force is greatest, approximately 
halfway through the contact interval. The kinetic energy of the impactor decreases rapidly 
as the impactor slows during contact with the target. Eventually, at tC/L = 0.25, the 
impactor velocity (and therefore its kinetic energy) decreases to zero and the elastic 
rebound begins. The kinetic energy of the impactor never returns to its initial level 
because approximately 80% of the energy has been transferred to the target in the form of 
strain energy and kinetic energy. The strain and kinetic energies in the target both 
increase rapidly during the contact with the impactor and remain constant after contact 
has ended (tC/L > 0.4). Both strain and kinetic energies maintain equal and constant 
values until the compressive stress wave generated by the impact reaches the far end of the 
free— free bar (tC/L = 1.0). A tensile stress wave is generated when the compressive pulse 
reflects from the stress free boundary [5]. The superposition of the incident and reflected 
pulses momentarily leaves the bar stress— free which causes the strain energy to decrease to 
zero. The kinetic energy simultaneously increases, maintaining a conservation of total 
energy. The reflection process is repeated at tC/L = 2.0, when the reflected pulse returns 
to the other end of the bar. Similar energy dissipation diagrams may prove useful in 
analyzing dynamic failure of more complex structures. 


Composite Plate 

The low velocity transverse impact of a composite plate made from Scotchply 
1003 prepreg [6] is now analyzed. The problem is depicted schematically in figure 5, and is 
described in detail in references [7,8]. A modified Hertzian contact stiffness has been 
proposed [9] for application to composite materials. Specifically, equation (6) is replaced 
by 


kt = E 33t (18) 

where E 33t is the transverse modulus of the plate. Plate membrane and bending stiffness 
material properties were calculated using the COBSTRAN (Composite Blade Structural 
Analyzer) computer code [10] which calculates elastic moduli of composite materials from 
known constituent properties and laminate ply orientations. 

A uniform square mesh of QUAD4 elements was used to model the 15.24 cm * 
15.24 cm (6 in * 6 in) target plate. A mesh convergence study was performed to establish 
the degree of mesh refinement necessary to arrive at a numerically converged solution. 
Three different meshes were considered, 25 * 25, 49 * 49, and 61 * 61 elements. Of these, 
the latter two produced essentially the same strain response for a given impact velocity and 
were therefore considered to be converged solutions. The results presented here were 
therefore calculated using the 49 * 49 element model. Five degrees of freedom (u x , u y , u z , 
0 X and 0 y ) were used at each nodal point, giving the model a total of 11510 degrees of 
freedom. The problem was solved on a Cray XMP in 52 CPU minutes. 
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The impactor used in the tests [7.8] was a uniform 2.54 cm (1 in) long, 
blunt— ended steel rod of radius 0.047625 cm ( 3 /i6 in). In the analysis a contact radius of 
0.047625 cm ( 3 / i6 in) was assumed in the Hertzian contact stiffness calculations. The 
calculated impact force history is shown in figure 6. Although no direct measurement of 
the impact force was obtained experimentally, the contact time was measured [8] and found 
to be 204 microseconds. This is in good agreement with the calculated result. Figure 6 
also shows that a secondary impact occurs during the latter half of the contact interval (t 
= 175 //sec), probably due to the vibration of the target plate during contact with the 
projectile. 


The resulting displacement response of the plate is shown in figure 7, where it 
has been assumed that no damage occurs in the target during contact with the impactor. 
This assumption is valid based on the available test data. Ultrasonic C— scans of the 
specimens after impact indicate that this level of impactor kinetic energy (10 Joules) is 
very near the threshold energy level required to cause damage [8] in specimens of this 
layup. As a result, very little damage occurs at this impactor velocity. The anisotropic 
bending stiffness of the target (figure 5) is evident from the elliptical displacement 
contours, as the flexural disturbance travels faster in the stiffer direction (figure 7). 

The strain response at gage A is compared to the calculated response in figure 8. 
The two curves are similar in amplitude and duration but the calculated strain appears to 
lag the measured values by approximately 25 microseconds. This may be due to the 
difficulty in establishing experimentally the precise time at which contact occurs based on 
strain gage readings taken at some distance from the point of contact. The comparison 
shown in figure 9 for gage B likewise shows a time shift of approximately 25 microseconds 
between the measured and the calculated response. The amplitude and duration of the 
calculated strain response correlate quite well with the measured signal. 


SUMMARY 

A simple means of modeling low velocity, non— damaging impact using 
NASTRAN was demonstrated. A nonlinear elastic contact model was included in the finite 
element analysis using NASTRAN transfer function definitions and nonlinear analysis 
capabilities. The same contact law was used to define the force— indentation relationship 
for two different impactor/target combinations. Results in both cases showed that the 
impact force and resulting transient structural response of the target compared well with 
experimentally measured values. 
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Figure 2: Impact Force for Bar Problem 
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Figure 5: Composite Plate Impact 
Specimen Configuration 
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Figure 7: Calculated Displacement Response for Composite Plate 
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APPENDIX 
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DAN TROWBRIDGE 


ANALEX - STRUCTRAL MECHANICS BRANCH 


ID TRANS. LOAD 
APP DISP 
TIME 60 
SOL 9 
CEND 

TITLE - COSMIC: TRANSIENT RESPONSE ANALYSIS: HERTZIAN IMPACT FF 
SUBTITLE - 36" AL. ROO 5/8 STEEL BALL V0-62.1 IN/S 

LABEL - ROO ' . < IMPACT 

* NONLINEAR LOAD 

NONLINEAR - 5 

$ INITIAL CONDITIONS SET 

IC - 1 

TFL-111 

SPC - 4 

TSTEP - 7 

$ OUTPUT STUFF 

SET 30 - 1.72,73.999,1001 
NLLOAD - 30 
STRESS (PR I NT) - 30 
DISP(PRINT) - 30 
BEGIN BULK 

$ 

$ EXTRA POINT m INDENTATION 

EPOINT, 1001 

GRID, 999. .-0.3125.0.0.0.0 
GRID. 1 . .0.0. 0.0. 0.0 

-O*«),.(1),-..(0.25).— 

CROO . 1 . 1 . 1 . 2 
-(143).»(1) ,»,»(1 ) . «(1 ) 

$ LUMP MASS OF IMP ACTOR 

CONM2, 200. 999, 0.9. 587-5, 0.0. 0.0, 0.0. , 4CON2-2 

+CON2-2.3.745-6, .3.745-6. . .3.745-6 

$ MATERIAL PROPERTIES 

PROO. 1.1 1.0. 196, 6. 14-3. 0.25 

MAT1 .11 .10.0+6, .0.33,2.5-4, . . ,+MATI-l 

+MAT1— 1 .35.0E6.36 . 0E6 , 27 . 0E6 

$ BOUNDRY CONDITIONS 

SPC 1 . 4 . 23456 , 1 . THRU . 1 45 

$ REMOVE DEGREES OF FREEDOM FROM IMPACTOR 

SPC1 , 4. 23456, 999 

$ TRANSFER FUNCTION TO DEFINE INDENTATION 

TF. 111. 1001. 0.+1. 0,0. 0,0.0., ,+TF— 1 

+TF-1 .999.1.-1 .0.0. 0,0.0, . . ,+TF— 2 

+TF— 2 .1.1. 1.0, 0.0, 0.0 

$ TIMING 

TSTEP . 7 . 2500 . 2 . 0-7 . 25 

S LOAD DEPENDENT ON DISPLACEMENT OF IMPACTOR 

N0LIN3. 5. 1. 1. 6.24+6. 1001. 1. 1.5 
$ SLOW DOWN IMPACTOR 

N0LIN3. 5. 999, 1. -6.24+6. 1001, 1. 1.5 

$ INITIAL CONDITIONS: IMPACTOR VELOCITY - 62.1 IN/SEC 

TIC. 1 .999.1 .0.0.62.1 

ENDDATA 
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ANALEX - STRUCTRAL MECHANICS BRANCH 


ID IMPACT. PLATE 
APP DISP 
TIME 120 
SOL 27 
CENO 

TITLE - IMPACT OF PLATE 49X49 : CENTERED ELEMENT 
SUBTITLE - TRANSIENT ANALYSIS: FIXED-FIXED: NO SYMMETRY 
SPC - 1 
IC - 3 

NONLINEAR - 5 
TSTEP - 1 
TFL - 111 

SET 15 - 999.2525.2526.2625.2626 
NLLOAD - 15 

SET 20 - 2525,3325.4125 
STRESS - 20 
BEGIN BULK 

$ • ••• EXTRA POINT TO HOLD INDENTATION •••*•• • •• 

EPOINT, 10001 

$ .... IMPACTOR ••• 3/8 IN DIAMETER • •••••••• 

GRID. 999, ,0.0.0. 0.-0. 1875 

CONM2. 200, 999. 0.8. 096-5, 0.0, 0.0. 0.0, , +CON2-2 
+CON2-2 , 7 . 459-6 . . 7 . 459-6 . . . 1 . 423-6 


$ ...... GRIDS AND CQUAD4 ELEMENTS DEFINING THE PLATE GO HERE ... 

$ MATERIAL PROPERT I ES . . . MAT2 CARDS GENERATED BY COBSTRAN 


PSHELL.1 .101 .0.15.201,1.0 

MAT 2 ,101 , 4 . 3E+06 , 2 . 9E+05 ,— 1 . 7E-03 , 2 . 8E+06 .-3 . 4E-02 . 5 . 7E+05 , 1 . 8E-04 ,+A101 
+A101 , 5 . 8E— 06 , 8 . 9E-06 , 5 . 0E— 1 3 

MAT 2 .201.5. 7E+06 , 2 . 9E+05 . -1 . 9E-04 . 1 . 4E+06 . -3 . 8E-03 . 5 . 7E+05 
$ BOUNDRY CONDITIONS 


SPC1 , 

1. 123456, 

101 . 

THRU. 150 

SPC1 . 

1, 123456, 

5001 

. THRU, 505e 

SPC1 . 

1. 123456. 

101 


-48* 

• 100 



SPC1 . 

1. 123456, 

150 



• 100 



-48 
SPC1 , 

1, 12456, 

999 



GRDSET 6 

$ TIME STEP INFO 

TSTEP. 1.2000. 1.0-7, 10 

$ LOAD DEPENDENT ON RELATIVE DISPLACEMENT OF IMPACTOR 


NOLIN3. 5,2525. 

3.+1 .945+5 

, 10001, 

0. 

1 .5 

NOLIN3. 5,2526, 

3.+1 .945+5 

, 10001, 

0. 

1 .5 

NOLIN3. 5.2625. 

3. +1 ,945+5 

, 10001. 

0, 

1 .5 

NOLIN3. 5.2626, 

3. + 1 .945+5 

. 10001. 

0. 

1.5 

$ 

SLOW DOWN 

IMPACTOR 



NOLIN3. 5. 999, 

3.-7.779+5 

. 10001. 

0. 

1 .5 


% TRANSFER FUNCTION TO CALCULATE INDENTATION 

TF, 111 .10001 ,0,+1 .0.00.0,00.0, . ,+TF-l 

+TF-1 ,999.3.-1 .0.00.0,00.0. . . .+TF-2 

+TF-2. 2525. 3. +0.25, 00. 0.00.0, . . .+TF-3 

+TF-3. 2526.3.40. 25.00.0,00.0, , , .4-TF-4 

+TF-4. 2625. 3. +0.25. 00. 0,00.0, . . ,+TF-5 

+TF-5. 2626. 3. +0.25, 00. 0.00.0 

$ INITIAL CONDITIONS: IMPACTOR VELOCITY - 1470 IN/SEC (122.5 FT/SEC) 

TIC. 3. 999, 3. 0.0, 1470.0 

ENDDATA 
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